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ABSTRACT 

A collisionless shock may be strongly modified by the presence of neutral atoms 
through the processes of charge exchange between ions and neutrals and ionization of 
the latter. These two processes lead to exchange of energy and momentum between 
charged and neutral particles both upstream and downstream of the shock. In par- 
ticular, neutrals that suffer a charge exchange downstream with shock-heated ions 
generate high velocity neutrals that have a finite probability of returning upstream. 
These neutrals might then deposit heat in the upstream plasma through ionization 
and charge exchange, thereby reducing the fluid Mach number. A consequence of this 
phenomenon, that we refer to as the neutral return flux, is a reduction of the shock com- 
pression factor and the formation of a shock precursor upstream. The scale length of 
the precursor is determined by the ionization and charge exchange interaction lengths 
of fast neutrals moving towards upstream infinity. In the case of a shock propagating 
in the interstellar medium, the effects of ion-neutral interactions are especially im- 
portant for shock velocities < 3000 km s^^. Such propagation velocities are common 
among shocks associated with supernova remnants, the primary candidate sources for 
the acceleration of Galactic cosmic rays. We then investigate the effects of the return 
flux of neutrals on the spectrum of test-particles accelerated at the shock. We find 
that, for shocks slower than 3000 km s~^, the particle energy spectrum steepens 
appreciably with respect to the naive expectation for a strong shock, namely oc E~^. 

Key words: acceleration of particles - atomic processes - line:profiles - ISM: su- 
pernova remnants 



1 INTRODUCTION 

The Physics of collisionless shock fronts is of paramount importance for numerous astrophysical sources, from supernova 
remnants (SNRs) and shocks in the Solar System to gamma-ray bursts and active galactic nuclei, as well as for the description 
of the process of cosmic ray acceleration. Such shocks are formed due to the mediation of electromagnetic instabilities while 
particle collisions are insignificant, thereby the name of collisionless. The shock itself, or the subshock in the case of a cosmic- 
ray (CR) modified shock, forms on a spatial scale of the order of ~ £,rL.th, where ^ > 1 and rL.th ~ W^'^ B^T^^ cm is the 
gyration radius of thermal particles at a temperature T — lO^Tg K and in a background magnetic field -B^ /iGQ The shock 
formation is associated with the development of electromagnetic streaming instabilities (e.g. Weibel instability), which can 
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^ Rigorously this statement applies to quasi-perpendicular shocks, while parallel shocks require a more detailed discussion. However for 
particle accelerating shocks, substantial magnetic field amplification must take place upstream, so that even if the shock configuration is 
initially parallel, due to compression, the magnetic field in the downstream plasma becomes mostly oblique. 
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amplify or even create magnetic fields (|Spitkovskvl2 008h. By definition, coUisionless shocks may only develop in ionized media. 
On the other hand, most astrophysical plasmas contain some fraction of neutral material, which will behave differently from 
the ionized component at the crossing of the shock. 

Ions are basically isotropized at the shock and heated to a temperature that mirrors the ram pressure of the incoming ion- 
ized fluid. If electrons are subject to the same fate, their temperature immediately behind the shock is m^/mp times the temper- 
ature of protons (ions). The temperature of electrons and pro tons may in fact turn out to be c loser to each other due to the elec- 
trom agnetic coupling between the two components (see e.g. Cargill fc PapadopoulosI 19881 : Ghavamian. Laming fc Rakowski 



20071 , and references therein). In addition, Coulomb scattering will also act towards temperature equilibration. This last pro- 



cess however, typically acts on much longer time-scales, which in the case of several observed SNRs (see e.g. lEllison et al, 
actually exceed the age of the source. 



As to the neutrals, these are insensitive to the shock transition and, to zeroth order approximation, they can cross the 
shock surface without interacting with the ions. The interaction of neutrals and ions occurs mainly through charge exchange 
and ionization and can change the structure of the shock rather dramatically because of the energy and momentum deposition 
that is involved. 

The cross section of these processes is of order ~ 10"^^ cm^, therefore, even in relatively tenuous plasmas, the rate at 
which they occur can lead to phenomena of potential astrophysical importance. Both charge exchange and ionization rates 
are proportional to the relative velocity between neutrals and ions. 

Upstream of a classical shock one expects that neutrals and ions are in some sort of local thermal equilibrium (as in the 
unshocked ISM), therefore they will have the same velocity and the same temperature. In these conditions there is neither net 
ionization, nor net charge exchange: clearly charge exchange will occur anyway because of the different velocities of neutrals 
and ions in different parts of their respective thermal distributions, but without changing the particle distributions. On the 
other hand, the neutrals crossing the shock front experience a thermal bath of hot ions with a bulk motion which is slower than 
their own (~ 4 times slower for a strong shock, if not modified by CRs). Moreover, while neutrals remain at their upstream 
temperature, ions are heated in the shock transition. As a result, a net velocity difference arises between the two components 
and both ionization and charge exchange are turned on, namely become effective at modifying the distribution functions of 
both species. 

Therefore, when a shock propagates into a partly neutral medium, an interesting chain of processes develops. When fast, 
cold neutrals undergo charge-exchange interactions with the slow hot ions downstream of the shock, some fraction of the 
resulting hot neutrals can cross the shock and move towards upstream infinity: a "return flux" develops. The relative velocity 
between these hot neutrals and the upstream ions triggers the onset of charge-exchange interactions that lead to the heating 
and slowing down of the ionized component of the upstream fluid. The system tends to develop a shock precursor, in which 
the fluid velocity gradually decreases from its value at upstream inflnity. As soon as the ions develop a velocity gradient in 
the upstream, charge exchange interactions become effective at modifying the particle distribution functions here as well, in 
a complex non-linear chain, in which the information is carried from the downstream to the upstream by the return flux of 
neutrals, which acts in a similar way as to what happens for CRs in modified shock theory (though see the discussion in 
for a comparison between the nature of CR- and neutral-induced precursors). The consequent reduction of the shock Mach 
number has potential implications, in turn, on the physics of dissipation and particle acceleration. 

The problem of describing a shock transition occurring in a partially neutral fluid can be treated in different ways but 
the general solution is known to consist of a final state (downstream infinity) in which neutrals eventually disappear through 
ionization and the ion density increases. This final state can be determined by adopting a fiuid approach bet ween upstream 



infini ty and downstream infinity and simply assuming conservation of the fluxes of mass, momentum and energy (|Morlino et al 



It is however not possible to use a naive two (or multi) fluid approach to adequately describe the complex system of 



neutrals and ions. This is because neutrals can only be taken to behave as a fluid on scales that are usually much larger 
than those of phenomenological importance. On such large scales, the system behaves as if the shock had become coUisional 
(the collisions being represented by charge exchange interactions) : also the neutrals are then "shocked" , through a transition 
region that has now become much thicker. At the same time, on a comparable length-scale, the neutral component gradually 
disappears because of ionization. 

Studies on this subject traditionally focused on the solar wind termination shock, but it is easy to see how the problem 
is extremely relevant also in the context of supernova remnant shocks. These shocks sometimes propagate in the cold, weakly 
ionized ISM, as shown by the associated Ha emission. In addition, several remnants showing an interesting phenomenology 
are actually thought to be interacting with dense cold clouds of quasi-neutral material. 



The Ha emission provides a powerful diagnostic tool for the conditions at the shock. As first recognized bv' Chevali er fc Raymond 



( 19781 ) and observed by Chevalier et al. (1980), the Ha profile detected in association with SNR shocks usually consists of two 
components, a narrow one, whose width refiects the temperature of the upstream medium, and a broad one, due to neutrals that 
have undergone charge-exchange with the hot downstream protons. This second component represents a uniq ue tool to measur e 
the temperature of ions, u sually very difficult to access otherwise. After the pioneering work of Chevalier fc Ravmond ( 19781 ) 
and Chevalier et al. I l|l980l ) several authors have further refined the use of Balmer emission as a diagnostics for SNR shocks. 
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Among ot hers, Ghavamian et alj l|2001l ) incorporated Monte Carlo modeling of Lyman a absorption, w hile Heng fc McCrav 
2007[l and lHeng et alj ([2007) included a more careful evaluation of the charge-exchange cross section and Ivan Adelsberg et al 



20081 ) considered the effect of multiple charge exchanges on the distribution function of the hot neutrals. However, all these 



works present limitations in the calculation of the distribution function of neutrals and do not take into account the effect of 
the return flux. 

The final goal of our investigation, to be illustrated in a forthcoming paper, is to use Balmer lines as a diagnostic tool for 
the detection of efficient CR acceleration in SNRs: if particle acceleration is effective, a sizeable fraction of the ram pressure 
is transformed into non-thermal particles, thereby reducing the heating at the shock surface. This should reflect in a lower 
ion temperature downstream of the shock, and hence in a narrower broad Balmer line. On the other hand, efficient particle 
acceleration also leads to the formation of a precursor upstream of the shock, where charge exchange can then occur. This 
leads to ion heating upstream, and correspondingly to a broader narrow Balmer line. ; Wagner et al. ( 2009 ) were the first to 
incorporate ionization, emission, and heating in a CR precursor into the mo deling of the Ha emission treating both CRs 
and neutrals as fiuids. A similar effort has be done by iRavmond et al.l (|201ll ) assuming a parametric structure for the CR 
precursor, fnt erestingly, s everal signatures of phenomena induced by the CR acceleration have been collected in the last few 
years (see e.g |Heng||2009l . for a review), but the absence of an appropriate mathematic tool still prevents us from describing 
them in a physically satisfactory way. Providing such a tool is the purpose of our series of papers, the present one being the 
first, while other implementations, like the inclusion of dynamically relevant CRs, will be the scope of forthcoming papers. 

The paper is organized as follows: in §[2] we write down the mathematical formalism for the solution of Vlasov equation for 
the neutrals and of the fluid equations for the ions. In §[3] we illustrate the main results in terms of modified shock dynamics, 
test particle acceleration in the presence of neutrals and Balmer emission. We conclude in § |4] 



2 THE MATHEMATICAL APPROACH 



The distribution function /jv(u, x) of neutrals interacting with ions through charge exchange and ionization is described by 
the Vlasov equation 

dfN 



dt 



-f V ■ V/jv = I^Nfi — ft/iV, 



(1) 



where fi is the distribution function of ions, assumed here to be a drifting Maxwellian distribution at a temperature T 
dependent upon location, with a bulk velocity representing the local speed of the ion plasma (see Eg. 1151 below). In Eq. [T]we 
have introduced 



PN{i 



(2) 



(3) 



and Vrel ~ \v — wl. The quantity /3i represents the rate of charge exchange (cross section ace) and ionization (cross section 
Cion) of a neutral with velocity v, while /3jv is the rate of charge exchange of an ion that becomes a neutral. The cross sections 
are all functions of the modulus of the relative velocity between an ion and a neutral and are plotted in Fig. [T] The ionization 
cross section is shown as a dash-dotted line, aee is the charge exchange cross section for hydrogen atoms including all possible 
flnal states. In Fig. [T] the solid line refers to the cross sect ion when the final state is the ground state, n = 1, while the 
dashed line refers to charge exchange with final state n = 2 ( Barnett el al. 1990l ). Following Heng fc McCra^ ; (2007 ), in orde r 



, first proposed bv I Jane v fc SmithI (|l993l ) 



to obtain cross sections for n > 2 we use the scaling relation ace,n = (2/n) 
Hence the total charge exchange cross section can be written as: 

oc oo 
fee — 0-ee,n ~ <Jce,l + (Tee^ (2/n)^ . (4) 

n=l n=2 

For (Tion we consider here only the ionization due to collisions with protons. For the sake of simplicity, we neglect the 
contributions to ionization due to collisions with electrons or heavy ions. These contributions may actually be comparable to 
that of protons, and therefore relevant for the calculation of the total Balmer emission (see e.g. Figure 1 of iHeng fc McCrav 
20071 ). We will introduce them in future works. 

In solving Eq. [T]we assume that a stationary situation is reached, so that dfN/dt — 0. Moreover we restrict our attention 
to the case of a plasma moving in one direction, z, so that the problem reduces to describing the evolution in the directions 
parallel (||) and perpendicular (_L) to z. Eq. [l]is then reduced to: 

dfN 



(5) 



where we used the fact that there cannot be gradients of the distribution function in the directions perpendicular to the 
symmetry axis of the problem, z. 
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Figur e 1. Cross section for cliarge excliange with n = 1 liydrogen atoms (solid line) and with n = 2 hydrogen atoms (dashed 
line) tearnett el al.! '1990"). The dotted line is the cross section for ionization of neutral Hydrogen due to collisions with pro- 
tons |janev & Smith 1993). All cross sections used are provided by the International Atomic Energy Agency via the web-site 
http : / /www-emidis . iaea.org/ALADDIN/. 



We describe here a novel way to solve the Vlasov equation, based on decomposing the distribution function /at in the 



and clearly /]v = "^f. /)^'. Since the neutrals do not feel the shock directly, for them there is no distinction between upstream 



sum of the neutrals that have suffered 0, 1, 2, ... , k processes of charge exchange. Each distribution function is named 



and downstream except for the interaction with different populations of ions on the two sides of the shock front. Upstream 
ions and neutrals are assumed to start (at z — )■ — oo) with the same bulk velocity and temperature, therefore charge exchange 
occurs at equilibrium and the distributions do not change appreciably (but see discussion below). On the other hand, after 
entering the downstream, the neutrals experience a very different environment and both ionization and charge exchange occur 
effectively. 

Let us consider the equation describing /^', namely the distribution of particles that have not suffered any charge 
exchange. This has no source terms and is easily seen to be: 



(6) 



with solution: 



-oo, , ux) X exp 



(7) 



Quite obviously, the number density of neutrals that did not suffer any charge exchange decreases exponentially after one 
interaction length, as defined by the ratio between uy and the rate Pi. It is worth noting that, in a plasma at equilibrium, 
charge exchange does occur although the net number of neutrals in a given volume of phase space does not change. This 
means that the position z — — oo is somewhat arbitrary and simply needs to be chosen far enough from the shock surface. It 



also means that in order to obtain the actual solution of the problem, one has to sum up a sufficiently large number of 



(fc) 



(fc) 



In the absence of any net velocity difference between ions and neutrals the sum of all /^y 



ff^ (z — — OO, , ux) for any z. The equation describing any is: 



I dz 



must return identically 



(8) 



The formal solution of this equation, as derived through elementary methods is: 



dz (fe_ 
Pn 



\z)fi exp 



rlz" 



A exp 



ft (2 



(9) 



where the integration constant A has to be determined by using the boundary conditions of the problem. These are different 
for v\\ > and < 0. Let us first consider the case v\\ < 0: for z — > +oo (downstream infinity) the number of particles that 
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have sufTered k scatterings must vanish, therefore, since — J^^dz/v^^fii diverges, one has to require 



A = - 



'fiexp 



Vll 



(10) 



For ?;|| > 0, the integration constant can be determined by recaUing that at upstream infinity there cannot be particles that 
have already sufTered k scatterings and that move with > 0, which implies A = 0. 



It follows that the general solution of Eq. [8]for the partial functions is 



and 



■IN 



Ak) _ 

J N ^ 



dz' (k- 

vn 



fi exp 



dz" 



Pi 



<0, 



> 0. 



From Eq. [HI and for «|| = 0, one easily obtains 



L J z.v\ 



The global solution of the Vlasov equation can now be written as the sum of all the partial functions: 

oo 

/]v(2:,«||,«±) = ^/^^'(2:,'U||,'Ux). 



(11) 



(12) 



(13) 



(14) 



In general, however, a good approximation to the solution is obtained when a sufficient number of partial functions is taken 
into account. The needed number of partial functions is determined by the physical scales of the problem, as we discuss below. 

Let us now describe how all the relevant are computed and what kind of approximations enter the calculation. 

First, let us consider the ions. We assume that newly ionized particles reach local thermal equilibrium with other ions 
in a time very short comparable with all other time scales. This seems a reasonable assumption in a scenario in which 
electromagnetic interactions lead to thermal equilibration. The ion distribution is assumed to be a Maxwellian 



/i(2,«|| 



ni{z) 



■ exp 



(v\\ ~Vi(z)f +v'i 



Vth,i{zY 



(15) 



[tt Vth,i{z)'^f'^'^ 

where ni{z), vth,i{z) and Vi{z) are the density, thermal velocity and bulk velocity respecti vely of ions at the position z. This 
form of the ion distribution function allows us to calculate Pi analytically as described by Pauls et al. 1 19951 ): 

I3i{z,v^^,v±_) = ■mpni{z)at{U*)Ut, (16) 



where C/, 



'4, ,2 



+ (Vi — Vu 



+ and at — Ccc + ^i is the total cross section for charge exchange and ionization. We 



checked that this approximation leads to the correct result with an accuracy of order few percent. 

The assumption of rapid (local) thermalization of newly produced ions with the bulk of ions deserves some discussion: this 
assumption is based on the fact that the thermalization process proceeds through the excitation of electromagnetic instabilities 
which are usually rather fast compared to the other processes involved. A caveat is due, however. The fastest process is most 
likely that of isotropization of the velocities of the newly created ions, which would make their effective thermal velocity 
comparable with their initial bulk velocity. While downstream of the shock such effective temperature is not too far from 
the temperature of thermal ions, in the upstream plasma the newly born ions would end up having a temperature much 
higher than that of the cold ions in the ISM. Whether the two populations can indeed reach some sort of thermal equilibrium 
upstream within a convection time to the shock is an open issue, and a rather difficult one. We are not aware of any previous 
discussion of this problem in the existing literature Q therefore we adopt here the simplest assumption also made by all 
previous works on this topic. It is however worth keeping in mind that the heating in the upstream fiuid predicted by our 
calculations may be somewhat overestimated in case of partial rather than total local thermalization of ions. 

ik) 

As for the calculation of the coefficients /3]y > this is the most challenging part of the work from the point of view of 
computation time: these are multi-dimensional integrals to be calculated on a multi-dimensional grid of values of {z,v\\,vi_) 
and in general the functions are far from being Maxwellian distributions, as discussed in § [S] We compute the 
following an approximate but physically motivated procedure, which is discussed in detail in Appendix [XI 

We then need to describe the evolution of the ion component under the action of charge exchange and ionization. As 
stressed above, ions are assumed to behave as a fluid, therefore their dynamics is described by a set of conservation equations 
which read: 

d 



[piVi + F„. 



0, 



(17) 



^ The possibility of non- Maxwellian proton distributions downstream of the shock as due to charge exchange reactions was presented 
by [Raymond et al.l | [2008tl . 
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^ [p^vf + Pg,, + Fmom] = 0, (18) 



d_ 

dz 



-piVi H '^^—-Pg,iVi + Fen 

2 79-1 



0, (19) 



where Fmass = rrip J (fivv^^fN, Fmom = rnp J d^vv'^^fN and Fen = rnp J d^vv\\ (wy + v\)fN are the fluxes of mass, momentum 
and energy of neutrals along the z direction. Notice that the terms on the right hand sides of Eqs. 1171 1181 and [19] are zero 
because mass, momentum and energy of the whole system (ions plus neutrals) are conserved. It is also worth stressing that 
nowhere here we assume that neutrals behave as a fluid. This point is in fact crucial in order to obtain the correct solution of 
the problem: two and three fluid approaches fail to describe the Physics correctly on scales smaller than the charge exchange 
length, as can be easily understood given that the neutrals feel the presence of the shock only indirectly, through their 
interactions with ions. In particular, it can be shown that fluid approaches in a stationary situation leads in general to the 
formation of a shock in the downstream neutral fluid: this is not what happens in nature as confirmed by our calculations. 

In practical terms our calculation is carried out by following an iterative procedure: we start with neutrals and ions in 
thermal equilibrium at "upstream infinity" (which translates, computationally, into many interaction lengths upstream of the 
shock). The relative density of the two components is parametrized through the ionization fraction. We start by fixing the 
profile of density, pressure and velocity of the ions, so that an ordinary shock front is formed at z = as predicted by solving 
the Rankine-Hugoniot equations for ions alone. The corresponding are calculated by solving Eqs. (|lip and (|12|) above, so 
that a mass, momentum and energy flux of neutrals are determined. At this point, using the conservation equations, Eqs. (|17p . 
(|18|) and (|19|) . we derive an updated profile of the quantities describing the ions' dynamics, and the cycle is restarted. 

Before embarking in the detailed explanation of the results of the calculations, let us illustrate the basic Physics that is 
expected, starting from the downstream region, where it is simplest. We will then proceed to introduce the novel phenomenon 
of return flux which is of the highest importance for astrophysical applications. 

The cold, fast neutrals from the upstream region penetrate the downstream region where ions have been slowed down 
by the shock and correspondingly heated up. In this situation the velocity difference between ions and neutrals ignites the 
processes of charge exchange and ionization (they also occur upstream, see below). A fast neutral that becomes an ion through 
charge exchange and/or ionization delivers momentum and energy to the ions. In the absence of ionization, one would expect 
that a large number of charge exchange events leads to thermalization between ions and neutrals. On scales much larger 
than the interaction lengths for the relevant processes we can think of using a "black box" approximation, namely of writing 
the mass, momentum and energy conservation equations as between upstream and downstream infinity (subscript 1 and 2 
respectively). We would then find that the final state (at downstream infinity) corresponds to a total compression factor Rtot 
that is the one appropriate for a system with total density equal to the sum of the densities of charged and neutral atoms (pi 
and pN respectively). In the presence of ionization only charged particles will be left at downstream infinity, with a density 

Pi, 2 = Rtot {pi,l + PN,i). 

As we already mentioned, the heating of neutrals through charge exchange is a phenomenon of crucial importance in that 
it affects the width of Balmer lines emitted by these hot neutrals. The line width is, in turn, a direct measurement of the 
temperature of ions downstream of the shock, at least on a spatial scale which is comparable with the excitation one (which 
is also comparable with the ionization scale). 

Neutrals that suffer one or few charge exchanges downstream have a very anisotropic distribution function. Charge 
exchange with ions in the part of the distribution function with ny < generates new neutrals which move against the stream 
of incoming particles in the shock frame. 

A few considerations are useful at this point: for a strong shock, if not modified by CRs, the thermal distribution of ions 
is centered around a bulk velocity which is ~ Vah/4: and has a spread of the order of ~ {3/8y^^Vsh- It is therefore easy to 
visualize that in the shock frame there are many ions that have a negative speed in the z direction, namely that move against 
the stream. As long as these particles are charged this does not represent a problem, in that their gyration radius is very small 
and they remain behind the shock[f|. However, following a charge exchange event, an ion with < becomes a neutral that 
keeps moving in the same direction. 

These neutrals have a finite probability of reaching upstream because again they are insensitive to the electromagnetic 
fields at the shock. When one such neutral happens to cross the shock and undergo a new interaction in the upstream, the 
associated deposition of energy and momentum will cause the heating of the upstream fiuid. This is what we will refer to as 
the neutral return flux, a phenomenon that could not be found in previous calculations, either because carried out in the fluid 
approximation or because concentrated on the downstream part of the plasma. Energy and momentum deposition upstream 
occur on a spatial scale which is the minimum between the charge exchange and the ionization scales. Whether one or the 
other dominates depends mainly upon the shock velocity. The relative velocity between the returning neutral and an ion is 



In fact, even a small fraction of ions might recross the shock, as discussed bv lEdmiston. Kennel fc Eichled l ll982l) . and eventually trigger 
the beginning of the injection process. 
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Figure 2. Density of particles in the individual distributions fj^ with k = 0, 1, 2, ... as indicated. The thick solid line is the density of 
neutrals as a function of the distance from the shock upstream (z < 0) and downstream (z > 0). 

of order ~ 2Vsh- For shocks with Vsh ~ 3000 km s~^ the relative velocity is such that the cross section for charge exchange 
is suppressed and the returning neutrals heat up the upstream fluid on the ionization scale (see Fig. [1} . For slower shocks 
the energy deposition takes place on the scale of charge exchange, since ionization is a threshold process in terms of relative 
speed. However one should also keep in mind that for high velocity shocks for most neutrals the first interaction downstream 
is of ionization type, therefore in this case the return flux is suppressed. 

We stress that, since the upstream gas is cold (r~ 10'' K), even a relatively small energy deposition may lead to a large 
increase in the ion temperature and a corresponding reduction in the Mach number of the plasma immediately upstream of 
the shock, which is thereby weakened. In ij 13. 21 we will describe how the return flux may affect the spectra of CRs accelerated 
at such type of shocks. 

3 RESULTS 

In this section we illustrate the main results of our kinetic calculations for a benchmark case with shock speed Vs — 2000 km s^^ , 
total density of 0.1 cm^'^ and 50 % ionization fraction (we specify the values of these parameters whenever they differ from 
the above). The temperature of the gas at downstream infinity is always assumed to be T = 10* K in order to be compatible 
with the presence of neutrals. The temperature of neutrals is assumed to be equal to that of ions, so at upstream infinity the 
two components are in thermal equilibrium. 

3.1 Dynamical properties of neutrals and ions 

In Fig. [2] we show the density of particles in each namely J d^vf^\v), as a function of the distance from the shock, 

both upstream (z < 0) and downstream (z > 0). As expected, the density of particles that did not suffer any charge exchange 
(k = 0) decreases monotonically with distance from upstream infinity. Clearly the concept of upstream infinity does not have 
physical relevance for us: moving it further away from the shock simply leads to requiring a larger number of to describe 
the total distribution of neutrals. Moving it closer to the shock surface is also possible, provided the neutrals that return from 
downstream after one charge exchange event become ionized or suffer another charge exchange within the distance associated 

(k) 

with the downstream infinity boundary. The density associated with the fj^ with > is not a monotonic function of the 
position because each distribution not only receives a contribution from the neutrals that have suffered (fc — 1) charge exchange 
reactions, but at the same time is also deprived of particles because of additional charge exchange and ionization interactions. 
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Figure 3. Left Panel: Distribution function of neutrals at the shock {z = 0, dashed) and at the locations zi = 2.1 X 10^^ cm (dotted 
line) and Z2 = 1.3 X 10^^ cm (dash-dotted line). The thick solid line is the distribution of ions immediately upstream of the shock. Right 
Panel: as in the left panel but for the downstream region. The thick solid line is the distribution of ions immediately downstream of the 
shock. 



One can see this phenomenon in Fig. [2] by comparing the contribution of fc = and fc = 1: the density of particles with k — 1 
increases with z within one interaction length of charge exchange, and then starts decreasing as a consequence of additional 
charge exchanges, which correspondingly contribute to the with fc > 1. Once the neutrals penetrate the downstream 
region they feel a larger density of ions and a larger velocity difference. The total density of neutrals (pjv) is plotted in Fig. [2] 
as a solid thick line. Approaching the shock from upstream, pN increases slightly because of the contribution of the return flux, 
namely of neutrals that come from the downstream region after a hot ion experienced a charge exchange thereby becoming 
a fast neutral that may move against the flow in the upstream. Moving from the shock vicinity to far upstream, the total 
density of neutrals becomes constant, because the returning neutrals disappear progressively, either due to ionization or to 
charge exchange. In the far downstream region the density of neutrals decreases again because of ionization. 

The distribution functions in phase space that result from our calculations are by no means Maxwellian functions, and 
contrary to previous calculations (e.g. see lHeng fc McCravl |2007)) they are computed at any point in space rather than being 
volume averaged. It is not easy to summarize in a plot the global properties of the functions /m since they depend on the 
position z and on velocity v — {vn,v±). In Fig. [3] we show /]v(z, , fx = 0) at two locations zi — ±2.1 x 10^^ cm and 
Z2 = ±1.3 X 10^" cm (the minus and plus signs refer to upstream and downstream respectively). The left and right panels 
show our results for the upstream and downstream section respectively. 

The thin and thick solid lines (both Maxwellians) in the left panel are the ion distribution at upstream inflnity and 
immediately upstream of the shock respectively. The long-dashed line is the distribution of neutrals at the shock location 
(since neutrals do not feel the shock directly their distribution is the same across the shock front). The short-dashed and dash- 
dotted lines refer to the distribution of neutrals at locations —zi and —Z2. At large distances from the shock the distribution 
of neutrals is identical to that of ions with temperature 10* K. In previous calculations of the structure of shocks in partially 
ionized media this situation was assumed to extend to the shock itself, namely nothing was happening in the upstream plasma. 
We flnd here that this is not a good description of reality because of the neutrals return flux. This component is clearly visible 
in the left panel of Fig. [3l in the region v^^ < 0. The most important effect induced by this return flux of neutrals is the heating 
of upstream ions. Due to charge exchange interactions between these hot upstream ions and cold upstream neutrals (the ones 
with > 0), the distribution function of the latter gets broadened, compared to the narrow Maxwellian at upstream infinity 
(see short-dashed and dash-dotted lines in the left panel of Fig. [3]). Clearly far enough upstream the return flux (wy < 0) 
disappears as a result of ionization and additional charge exchange reactions, so that, as we stressed above, the distribution 
of neutrals at upstream inflnity sits on the ion distribution at T = 10* K. 

The situation downstream of the shock is somewhat simpler and is illustrated in the right panel of Fig. [S] The distribution 
of neutrals at the shock (dashed line) is the same as in the left panel, but the neutrals at zi and Z2 developed a very broad 
distribution, roughly with the same width of the ions' distribution immediately behind the shock (thick solid line). This is 
the result of efficient charge exchange between the cold (T = 10* K) neutrals and the hot ions behind the shock. Nonetheless, 
one can also see that at the distances zi and Z2 considered here there is still a leftover narrow distribution of particles 
from upstream. Their contribution to the total number and energy is however rather small in the downstream plasma. 
Moving towards downstream infinity the neutrals disappear as a result of ionization, although for slow moving plasmas this 
phenomenon can occur very far from the shock because of the small ionization cross section. 
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Figure 4. Le/i Panel: Density of neutrals (lower curves) and ions (upper curves) as a function of the distance from the shock upstream 
and downstream, for different values of the shock velocity. Right Panel: Normalized ion velocity for different values of the shock speed. 



The role of the return flux is better illustrated by the dynamical quantities associated with ions. In Fig. |3] we show the 
density (left) and velocity (right) of the ion plasma as a function of the position across the shock {z = 0). The different 
curves refer to different values of the shock velocity as labelled. One can see that the density stays constant upstream until 
very close to the shock front, where the return flux becomes important. In this region, the fast neutrals from downstream 
deposit energy into the ion plasma through both ionization and charge exchange, thereby heating the gas. These reactions also 
deposit momentum in the —z direction, thereby slowing down the ion plasma (right panel). The combination of this induced 
deceleration and of the ionization process causes the ion density to increase (left panel). One flnal comment is deserved by the 
right panel of Fig. [l] here we see that with decreasing shock velocity between 3000 and 1500 km s~^ the precursor becomes 
progressively more pronounced. This trend suddenly changes for Vsfi=1000 km s~^, when the precursor becomes shorter and 
steeper: this is due to the fact that in this case the return flux is destroyed by charge-exchange events rather than by ionization, 
and those occur on a shorter scale. 

It is important to realize that the whole system of neutrals and ions on very large scales, namely between upstream 
infinity and downstream infinity, must behave as a black box in which the standard conservation equations must apply (a 
generalized version of this statement can also be written in the presence of CRs). Since for T = 10^ K we have a sonic Mach 
number M = 85 Vsh/ {WOO km s~^) 3> 1, the shocks we are dealing with are indeed strong: therefore, the velocity of ions and 
neutrals, at upstream infinity, has to be ~ 4 times larger than the velocity of the far downstream plasma, made of ions only 
because of ionization. This fact is clearly visible in the right panel of Fig. [4] From mass conservation one immediately obtains 
that 

PiVsh + pNVsh = ^-^PF ^ Pf = 4 (pi + pn) , (20) 

where pp is the density in the far downstream plasma. For equal densities of neutrals and ions at upstream infinity, pp/pi = 8. 
This is clearly visible in the left panel of Fig.[4l where the upstream density is 0.05 cm"'^ and the density at downstream infinity 
pp is 0.4 cm^^, independently of the shock velocity. Increasing the shock velocity leads to reducing the relative importance 
of the return flux, as illustrated in Fig. [5] where we plot the temperature of ions as a function of location around the shock. 

All curves start at T = 10^ K at upstream infinity, but it is clear that, on the way to the shock, ions are heated to 
much higher temperatures due to the effect of the return fiux described above. The fact that the gas temperature is of order 
~ lO'^ K immediately upstream of the shock leads to a dramatic reduction in the ion Mach number and therefore in the 
compression factor at the shock which is appreciably smaller than the canonical value of 4, especially for low shock velocities 
(below ~ 3000 km s"^). 

In the perspective of discussing the implications of these physical processes for diffusive particle acceleration at the 
shock, it is worth stressing that the effect of neutrals on the shock structure is that of producing a precursor upstream of 
the shock, but that the nature of this p recursor is totally unrel ated to the CR induced precursor that is found in non-linear 
theories of particle acceleration (see e.g lMalkov fc Drurvll200ll . for a review). The neutrals-induced precursor is simply due 
to the fact that neutrals can carry information from downstream to upstream of the shock, and deposit such information on 
spatial scales which are fully determined by the cross sections of charge exchange and ionization. In the case of CR-induced 
precursors, the plasma upstream is slowed down by the pressure of accelerated particles, and the spatial scale of the precursor 
is determined by the diffusion properties of the medium and by the spectrum of accelerated particles (in general the precursor 
is larger whenever the spectra are harder, which reflect more efficient acceleration). All these properties are poorly known. 
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Figure 5. Temperature of ions as a function of position around the shock, for different values of the shock speed. 



while the spatial extent of the neutral-induced precursor is basically fixed by the cross sections for charge exchange and 
ionization. The heating of the plasma in a CR-induced precursor is due to adiabatic heating and to turbulent heating. The 
latter depends on unknown details of wave damping on the gas, and its efficiency can only be parametrized. The process 
cannot be too effective otherwise the wave amplification that is responsible for effective diffusion upstream is inhibited. In 
a neutral-induced precursor, the heating is due to energy and momentum deposition of ions produced in charge exchange 
and ionizations reactions of returning neutrals upstream. The only uncertainty here is due to the unknown rapidity of ion 
assimilation in the thermal plasma (see also discussion below). 



3.2 Acceleration of test particles in partially ionized media 

In the context of diffusive shock acceleration the spectrum of test particles accelerated at the shock is a power law N{E) oc E~'' 
with a slope 7 fully determined by the compression factor at the shock: 

7^^. (21) 
r — 1 

For a strong shock (sonic Mach number Af ^ 1), if not modified by CRs, the compression factor r — >■ 4 and the spectrum 
reaches its asymptotic shape N{E) ~ 

As discussed in il3. li the presence of neutrals induces the formation of a precursor upstream of the shock front. The ion 
temperature immediately upstream of the shock may become 2-3 order of magnitude larger than the temperature at upstream 
infinity, hence the Mach number at the shock is much reduced. The importance of this effect depends upon the shock velocity. 
In the left panel of Fig. |6]we show the temperature immediately before (Ti) and behind {T2) the shock as a function of the 
shock velocity. The dotted line illustrates the downstream ion temperature in the absence of neutrals. 

The temperature immediately behind the shock is basically the same with or without neutrals. This fact can be understood, 
in a qualitative way, by considering a "black box" description of the whole system, namely writing the conservation equations 
for mass, momentum and energy between upstream and downstream infinity, while ignoring the detailed physics in between. 
When this is done, the ion temperature at downstream infinity is fixed, independent of the presence of neutrals. If one further 
considers that charge exchange processes in the downstream have very little effect on the temperature of ions, the above result 
is readily interpreted. 

What is most impressive is that the temperature of ions upstream grows to very large values, Ti ~ 10® — 10^ K due to 
the presence of neutrals. We stress that such a heating may be much stronger than any other mechanism, like for instance 
the turbulent heatin g due to the damping of Alfven waves, which may be expected to be effective in CR-modified shocks (e.g. 
CaprioU et"Zll2009l . and references therein) . The reason is that the reservoir of energy to be damped by Alfven heating is the 



energy in the form of magnetic turbulence (typically less than 1 per cent of the bulk energy upstream), while here the neutral 
return flux itself accounts for a potentially large fraction of the bulk pressure (see the right panel in fig. |4]). 
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Figure 6. Left Panel: Temperature of ions immediately upstream (Ti) and downstream (T2) of tlie sliock. Tiie dotted line is the 
downstream temperature in the absence of neutrals. Right Panel: Mach number of the ions fluid at upstream infinity (solid line) and 
immediately upstream of the shock (dashed line). The dash-dotted line shows the compression factor at the shock in the presence of 
neutrals. 

Ti grows with shock velocity until it reaches a maximum around ~ 2000 km s~^, and decreases for faster shocks. This 
trend reflects the physical essence of the return flux: for shock velocities smaller than ~ (2 — 3) x 10^ km s~^ neutrals entering 
the downstream plasma are most likely to suffer a charge exchange interaction, and the resulting fast neutral has a finite 
probability of having D|| < thereby contributing to the return flux. As we already mentioned, in this situation the ion plasma 
upstream gets heated and a precursor is formed. On the other hand for V^h > (2 — 3) x 10^ km s~^, ionization occurs before 
charge exchange and the return fiux is correspondingly suppressed. This explains the decline of Ti in the left panel of Fig. [6] 
in the high shock speed region. 

In the right panel of Fig. [6] we plot the Mach number at upstream infinity (Mo, solid line), the Mach number of ions 
immediately upstream of the shock {Mi, dotted line) and the compression factor at the (sub)shock, Rsub (dashed line), as 
functions of the shock velocity. For all values of Mo the compression factor that would be derived in the absence of neutrals 
is ~ 4, but the action of neutrals is such that the compression factor drops to values below 2 for Vsh < 1500 km s~^ and 
gradually grows to 4, which is reached however only for Vsh ~ 3000 km s~^. 

Even on a qualitative basis it is clear that the presence of neutrals, by affecting the compression factor at the shock, will 
also affect the spectrum of test particles accelerated at shocks: more specifically, the latter will become steeper than standard 
test particle theory would predict for a high Mach number shock. It is also clear that the spectrum of accelerated particles 
must be concave to some extent because the compression factor experienced by low energy particles is closer to Rsub, while 
higher energy particles experience a compression factor closer to the standard one, ~ 4. 

In order to estimate this effect we introduce the energy-dependent compression factor 

R{E) = (22) 

where 

n- = n., + ^/rf.giV(i.,.) (23) 

represents an effective fluid velocity upstream (1) and downstream (2) of the shock, as experienced by particles with energy E, 
No{E) is the spectrum of particles at the shock location, and ui (112) is the fluid velocity immediately upstream (downstream) 
of the shock. The spatial integral is extended to upstream infinity for and to D{E)/u2 downstream for u^^\ where D{E) 
is the diffusion coefficient for a particle with energy E. The reason is that the downstream region in principle extends to 
infinity, but the particles that can return to the shock due to diffusion are only those that reside within a region of size 
D{E)/u2 downstream. All the particles in the upstream region are eventually advected towards the shock front. Here the 
diffusion coefficient is assumed to be Bohm-like with a magnetic field of ^ 10 fiG upstream and ~ \/TT x 10 /iG downstream 
(formally the factor y/Tl holds only for compression of a turbulent magnetic field at strong shocks with compression factor 4, 
but this is not very important in this context). 

The compression factor R{E) provides an estimate of the actual compression factor experienced by particles with energy 
E. The slope of the spectrum is therefore defined as ^{E) = (see Eq. I21|) and plotted in Fig. [7] for four values of the 

particle energy {E = 1, 100, 100, 1000 GeV) as a function of the shock speed. 

Very high energy particles {E — 1 TeV, dotted line) sample almost the entire large-scale structure of the shock so that 
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Figure 7. Slope of the spectrum of accelerated test particles for E = I, 10, 100, 1000 GeV, as a function of the shock speed. 



for them the effective compression factor is close to 4. The corresponding spectral slope varies between 2 (for very slow and 
very fast shocks) and 2.1 for Vsh ~ 1500 km s~^. 

For particles with E — 1 GeV, the slope is considerably affected by the presence of neutrals, becoming as large as ~ 4.5 
for Vsh ~ 500 km s~^. The slope approaches the canonical value of 2 only for Vsh ~ 4000 km s"'^. The effect of neutrals is 
very evident also for 100 GeV particles (short-dashed line): the slope gets as steep as ~ 2.7 for Vsh ~ 1000 km s~^, and is 
always larger than 2.3 for Vsh ~ 2500 km s^^. 

These results clearly show how the spectrum of accelerated particles is affected in a very important way by the presence 
of neutrals for ionization fraction of 50% (our benchmark case) and shock velocity J; 4000 km s~^. In Fig. [8] we plot the 
spectral slope of test particles for E = 1, 10, 100, 1000 GeV, Vsh = 2000 km s"\ n = 0.1 cm ^ as a function of the fraction 
of neutrals. 

A departure of the spectral slope of accelerated particles from the canonical value of 2 is observed as soon as the neutral 
fraction is non-vanishing. The spectrum becomes especially steep at low energies, since these particles probe spatial scales 
that are entirely contained within the precursor induced by the return flux of neutrals rather than the global extent of the 
system. For a neutral fraction ~ 0.8 even the spectral slope at ~ 1 TeV is ~ 2.3. 

It is worth stressing the comparison between the spectral steepening induced by the presence of neutrals and that 
induced by non-linear effects in particle acceleration. Very efficient acceleration may lead to steep spectra at energies below a 
few GeV, as a consequence of the formation of a pronounced CR-induced precursor: the steepening is caused by the fact that 
low energy particles can feel the compression factor at the subshock, which may be weakened if acceleration is efficient. These 
effects disappear above a few GeV because h igher energy particles feel the whole precursor and their spectrum is appreciably 
harder than E~^ . It has been pointed out bv lCaprioli et al. I l|2O10l ) that spectra steeper than E ^ (and correspondingly lower 
efficiencies of particle acceleration) can be obtained even in the context of the non-linear theory of diffusive shock acceleration 
if the velocity of the scattering centers is taken into account. However it is worth recalling that in these cases the results are 
strongly dependent upon the detailed nature of the waves and on their helicity: in principle the same effect may lead to harder 
spectra rather than to a steepening. On the other hand, the neutral return flux induces a precursor whose length scale (the 
charge-exchange/ionization mean free path) is typically much larger than the diffusion length of few GeV particles, thereby 
potentially affecting several decades of the CR spectrum up to multi-TeV energies, as illustrated in Fig. |8l 

The possibility of obtaining particle spectra significantly steeper than E~^ well above 1 GeV due to the presence of a 
relevant neutral component is particularly intriguing because of the recent detection of many SNRs by Fermi and AGILE in 
the GeV band, and by HESS, VERITAS and MAGIC in the TeV band. These observations are showing compelling evidence 
of gamma-ray spectra typically in the range 

£;-2.2 _ ^-2.4 j-^j. shell-like SNRs (with the exception of RX J1713. 7- 3946 and 
Vela Jr.) and even steeper (_E~^'^ — E~^'^) for SNRs interacting with (partially neutral) molecular clouds (see e.g. ICaprioli 
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20111 . for a review and a wider discussion). Although challenging from the theoretical point of view, the presence of neutrals 
may play an important role in making this issue milder. 



4 CONCLUSIONS 

The structure of a coUisionless shock wave is profoundly affected by the presence of neutral atoms in the medium in which the 
shock propagates. The coupling between the shocked ions and the neutrals occurs through the processes of charge exchange 
and ionization and leads to strong modifications of the shock structure. 

Our calculations are based on a novel procedure that allows us to solve semi-analytically the Vlasov equation describing the 
behaviour of neutrals and the fluid equations describing the ionized plasma. Both ions and neutrals are evolved from upstream 
infinity to downstream infinity and the structure of the shock is calculated iteratively. The main physical phenomenon that 
we discovered, and the one that is responsible for all the interesting effects discussed here, is the existence of a return flux 
of neutrals in the upstream: a cold neutral crossing the shock from upstream may suffer a charge exchange reaction with a 
hot ion in the downstream; some of these interactions will involve ions that are moving towards the shock (negative value 
of the velocity component parallel to the shock normal, v\\); when this happens, the ion is transformed into a neutral that 
keeps moving with the same velocity, v\\ < 0, thereby recrossing the shock and reaching upstream. The return flux of neutrals 
created in this way is then dissipated upstream through additional charge exchange and ionization reactions. This leads to 
heating of the upstream ions and to the consequent decrease of the shock Mach number with respect to the value at upstream 
infinity. We find that in some cases the Mach number drops to values of order ~ 2 immediately before the shock. This implies 
that: 1) a precursor is induced upstream of the shock by the neutral return flux; 2) the compression factor at the shock is 
lowered much below the standard value of 4 that applies to strong shocks, if not modifled by CRs. 

There is an intrinsic velocity scale in the problem, which is of order 2000 — 3000 km s~^ established by the cross sections 
of ionization and charge exchange: for shocks with Vsh < (2 — 3) x 10^ km s~^, a neutral that crosses the shock is more 
likely to suffer a charge exchange reaction with an ion downstream rather than being ionized. In these conditions a flux of 
neutrals with < is created and the shock is profoundly modifled. At higher shock velocities the neutral gets ionized before 
it suffers charge exchange, therefore the return flux is suppressed and the shock structure is not affected appreciably. The 
spatial scale of the precursor induced by the return flux is determined by the shortest between the ionization and the charge 
exchange interaction lengths. This scale is again a function of shock velocity: the relative velocity between an incoming ion 
and a returning neutral is of order ~ 2Vsh- If ^Vsh < (2 — 3) x 10^ km s~^ the main process for dissipation of the return flux 
upstream is charge exchange, otherwise energy and momentum are deposited through ionization. 

In our calculations neutrals are described through the Vlasov equation, therefore we can determine their distribution 
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function at any location upstream and downstream. These distributions are not Maxwellian in shape: upstream of the shock 
the particle distribution function is bimodal, with a roughly Maxwellian peak that describes the neutrals that did not suffer 
charge exchange reactions yet, and a broader component in the region un < that disappears while moving towards upstream 
infinity. In the downstream region the distribution function competing to neutrals that have not suffered any charge exchange 
rapidly vanishes and is replaced by a broad component, made of hot neutrals and reflecting the charge exchange reactions 
with hot ions. Moving towards downstream infinity all neutrals eventually disappear because of ionization. The temperature 
of ions upstream is a strong function of distance from the shock, as a consequence of the return flux. Downstream the ion 
temperature varies little. 

The implications of the formation of a neutral-induced precursor for particle acceleration at collisionless shocks are 
of considerable interest. The spectrum of test particles resulting from diffusive shock acceleration is determined by the 
compression factor the particles experience while they diffuse in the region surrounding the shock. We find that in the 
presence of neutrals this compression ratio, even in the case of a high Mach number shock, is smaller than the canonical 
value of 4. It follows that for all particles for which the diffusion length upstream, ~ D{E)/Vsh is shorter than the precursor 
length, the spectrum is steeper than _E~^. Our calculations provide a quantitative confirmation of this qualitative expectation: 
assuming conditions typical of supernova remnant shocks, we obtain that the spectrum in the 1 — 10 GeV range may become 
as steep as ~ E''^ for V^h = 1000 km s"^; at 100 GeV the spectrum is still ~ E~^-^ and flattens to ~ E~'^'^ for TeV 
energies. The precursor weakens at larger shock velocities and the particle spectrum becomes progressively less deviant: for 
Vsh ~ 2000 km s~^, the spectrum is always between and E~^'^ . At even larger shock velocities the standard results are 

reproduced. 

These results are obtained by assuming that the accelerated particles behave as test particles. Their dynamical reaction 
may become important in realistic situations and this leads to a bunch of new efi'ects that will be discussed in a forthcoming 
paper. Most notably the accelerated particles induce the formation of a precursor upstream of the shock which slows down ions 
with respect to neutrals, thereby establishing a velocity difference that triggers charge exchange far upstream of the shock. 
This effect greatly modifies the shock structure as we discuss in a forthcoming paper, where the calculations illustrated here 
will be generalized to the case of non-linear shock acceleration. There, we will also illustrate the results of our calculations on 
the width of Bahner lines in shocks where efficient CR acceleration takes place. 
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APPENDIX A: CALCULATION OF (3 

(k) 

The calculation of the scattering rates has to be carried out in an approximate way since the calculation of the integrals 

(k) 



in d V on the entire three dimensional grid {z, t;|| 



leads to exceedingly long computation times. The expressions for / 



provided in Eqs. ((TTJ and (|12|) can be easily rewritten in an easy form separating charge exchange interactions occurring in 
the upstream and in the downstream. For instance, let us consider a point upstream of the shock, at 2: < 0, and let us consider 
first the case U|| < 0. Then, from Eq. (|ll|l one has: 
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where and f'^'^ are the contributions to f'^' deriving from charge exchange events occurred upstream (u) and downstream 
(d). Similarly, for > and 2 < 0, from Ea.[T2l 
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In Eq. (|A1|) one should notice that the fraction of ions in the distribution function fi with < is tiny, so that to a good 



approximation / 



< 0,2 < 0) 



J N,d- 



In other words the main contribution to the with un < upstream comes from 



charge exchange events that have occurred downstream. On the other hand, /^'(un > 0, 2 < 0) 



Ak) 
J N,m 



namely the part of 



the distribution function with iiy > is completely determined by the charge exchange events that occur upstream. Hence, the 
separation of into the two contributions /jv,u and /jv,ti roughly coincides with the separation of the distribution function 
into f^\z < 0,v\\ < 0) and fl^\z < 0,v\\ > 0). A similar line of thought in the downstream region leads to: 



/f,'>(2>0,^;|| <0) = /^ 



(k) 



(A3) 
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and 

/(,^)(^>0,^;|| >0) = (A4) 

As stressed in the main part of the paper and illustrated in Fig. [3l the distribution functions are not Maxwellians, neither 
are Maxwellian the individual fl^\ although for fc ^ 1 their shape eventually becomes closer to that of a Maxwellian. 
The functions and /^^^, however, are much more similar to a Maxwellian than the total distribution. Hence, in order to 

( k) (k) 

simplify the calculation of the scattering rates we determine the moments of fj^ ^ and |j and we determine their contribution 
to as the sum of the contributions of Majcwellians with the same moments as f^\^ and /^^. In other words, for f^^^ and 
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and we write: 
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where C/*,iv,u = \l ^I'l^Jjv.n'' + ("^^,1 " ^'ll)^ + «i and U,,N,d = J ^v[h,N,d + («^*d " + 
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